#install.packages("plantecophys")
library(plantecophys)
library(bigleaf)

#path to save image output
wd<-"~/HeatStress/Figures/"

#constant model parameters as defined in Supplementary Table 1
alpha=0.3  #reflectance of human skin
rho=1.2 #density of air (kg/m3) 
cp=1010 #specific heat capacity of air (Joules/kg/°C)
pa = 101000 #Pa; normal air pressure at sea level
es = 0.97 #Human skin emissivity (Brewster, M. Quinn (1992))  
hc = 3.3 #W m-2 K-1; natural convection when wind speed<0.2 m s-1, from literature survey shown in Supplementary Fig.7
Ts = 36  #°C; skin temperature in warm conditions prior to core body temperature rises (Fiala et al. 2001; Cuddy et al. 2014; ISO Standard Parsons 1999)
fs = 0.8 #the fraction of total skin area effectively involved in radiant heat exchange (taken as a constant 0.8 from Steadman 1979, Fiala et al. 1999)


#--------test the model against PSU heat data-----------
#See detailed description of the PSU HEAT experiments in Wolf et al. 2022. Critical environmental limits for young, healthy adults (PSU HEAT Project). Journal of Applied Physiology 132, 327–333. https://doi.org/10.1152/japplphysiol.00737.2021
#In the first three scenarios, the air temperature T is fixed at about 36°C, 38°C, and 40°C, respectively, 
#whereas the vapor pressure is gradually increased until the core temperature of the subject starts to rise. 
#In the other three scenarios, the vapor pressure pv is fixed at around 2.7 kPa, 2.1 kPa, and 1.6 kPa, respectively, 
#whereas the air temperature is gradually increased until the subject’s core temperature starts to rise.

#1) ==========MinAct experiment==============
MH_m <- c(77.9, 82.6, 80.4, 81.5, 87.7, 87.4)
MH_sd <- c(15.7, 17.1, 8.5, 11.4, 10.7, 12.3)

Ta_m = c(36, 38, 40, 43.8, 47.3, 50.6)
Ta_sd = c(0, 0, 0, 2.0, 1.6, 2.0)
#vapor pressure in the air (Pa); same as e_a=RH*esat_air
Pa_m = c(29.6, 29.8, 27.7, 20, 16, 12)*133.322 #mmHg=133.322 pa
Pa_sd = c(2.2, 2.3, 2.4, 0, 0, 0)*133.322
RH_m = c(66.4, 60.5, 50.2, 29.2, 20.3, 12.7)
RH_sd = c(5.4, 5.1, 4.3, 2.8, 1.5, 1.5)

#net longwave radiation at skin surface
Ln_m = fs*(es*5.67*10^-8*((Ta_m+273.15)^4 - (Ts+273.15)^4))
Ln_sd1 = fs*(es*5.67*10^-8*((Ta_m+Ta_sd+273.15)^4 - (Ts+273.15)^4))
Ln_sd2 = fs*(es*5.67*10^-8*((Ta_m-Ta_sd+273.15)^4 - (Ts+273.15)^4))

lambda0=latent.heat.vaporization(Ts)
gamma0=psychrometric.constant(Ts, pa/1000)*1000
esat_Ts=Esat.slope(Ts, formula="Alduchov_1996")$Esat*1000

#body surface energy balance
G1_m = Ln_m + MH_m - hc*(Ts - Ta_m) - hc/gamma0*(esat_Ts - Pa_m)
G1_sd1 = Ln_sd1 + (MH_m+MH_sd) - hc*(Ts - (Ta_m+Ta_sd)) - hc/gamma0*(esat_Ts - (Pa_m+Pa_sd))
G1_sd2 = Ln_sd2 + (MH_m-MH_sd) - hc*(Ts - (Ta_m-Ta_sd)) - hc/gamma0*(esat_Ts - (Pa_m-Pa_sd))
G1_m;G1_sd1;G1_sd2


#2) ============LightAmb experiment==============
MH2_m <- c(128.5, 136, 133.5, 130.3, 134.4, 137.3)
MH2_sd <- c(17.5, 15.8, 13.2, 10.3, 14.5, 18.7)

Ta2_m = c(36, 38, 40, 38.9, 43.0, 44.5)
Ta2_sd = c(0, 0, 0, 1.0, 1.5, 2.4)
#vapor pressure in the air (Pa); same as e_a=RH*esat_air
Pa2_m = c(24.6, 22.6, 20.7, 20, 16, 12)*133.322 #mmHg=133.322 pa
Pa2_sd = c(1.7, 1.0, 1.8, 0, 0, 0)*133.322
RH2_m = c(55.5, 45.5, 37.3, 38.4, 24.9, 18.6)
RH2_sd = c(3.9, 2.7, 3.3, 1.8, 2.0, 3.3)

#net longwave radiation at skin surface
Ln2_m = fs*(es*5.67*10^-8*((Ta2_m+273.15)^4 - (Ts+273.15)^4))
Ln2_sd1 = fs*(es*5.67*10^-8*((Ta2_m+Ta2_sd+273.15)^4 - (Ts+273.15)^4))
Ln2_sd2 = fs*(es*5.67*10^-8*((Ta2_m-Ta2_sd+273.15)^4 - (Ts+273.15)^4))

#body surface energy balance
G2_m = Ln2_m + MH2_m - hc*(Ts - Ta2_m) - hc/gamma0*(esat_Ts - Pa2_m)
G2_sd1 = Ln2_sd1 + (MH2_m+MH2_sd) - hc*(Ts - (Ta2_m+Ta2_sd)) - hc/gamma0*(esat_Ts - (Pa2_m+Pa2_sd))
G2_sd2 = Ln2_sd2 + (MH2_m-MH2_sd) - hc*(Ts - (Ta2_m-Ta2_sd)) - hc/gamma0*(esat_Ts - (Pa2_m-Pa2_sd))
G2_m;G2_sd1;G2_sd2


#====G model: add isopleth of G=0
#given Ta of 30 to 60, determine vapor pressure Pa and RH for G=0
Ta0 <- seq(30,60,1)
Ln0 <- fs*(es*5.67*10^-8*((Ta0+273.15)^4 - (Ts+273.15)^4))

#for MinAct: 0 = Ln0 + mean(MH_m) - hc*(Ts - Ta0) - hc/gamma0*(esat_Ts - Pa1)
Pa1 <- esat_Ts - (Ln0 + mean(MH_m) - hc*(Ts - Ta0))*gamma0/hc
esat_air=Esat.slope(Ta0, formula="Alduchov_1996")$Esat*1000
RH1 <- Pa1/esat_air*100
LE_req <- Ln0 + mean(MH_m) - hc*(Ts - Ta0) #latent heat of free evaporation under natural convection
#LE_req well below the maximum sweat capacity of 500/400 W m-2 or 1250/1000 g h-1 for acclimated/non-acclimated adults


#for LightAmb: 0 = Ln0 + mean(MH2_m) - hc*(Ts - Ta0) - hc/gamma0*(esat_Ts - Pa2)
Pa2 <- esat_Ts - (Ln0 + mean(MH2_m) - hc*(Ts - Ta0))*gamma0/hc
RH2 <- Pa2/esat_air*100
LE_req2 <- Ln0 + mean(MH2_m) - hc*(Ts - Ta0) 
#LE_req2 well below the maximum sweat capacity of 500/400 W m-2 or 1250/1000 g h-1 for acclimated/non-acclimated adults


#====add contour of TW
Pa0 <- seq(0, 6, 0.2)
Tw <- matrix(nrow=length(Ta0),ncol=length(Pa0))
for (n in 1:length(Ta0)) {
  for (m in 1:length(Pa0)) {
    esat_air0=Esat.slope(Ta0[n], formula="Alduchov_1996")$Esat
    VPD <- esat_air0 - Pa0[m] #kPa
    Tw[n,m] <- wetbulb.temp(Ta0[n], pa/1000, VPD, Esat.formula="Alduchov_1996") #degC
}}



#===main Figure 2
pdf(paste0(wd,"/Fig.2.pdf"), width = 9/2.54, height = 7.5/2.54 )

# windowsFonts(A = windowsFont("Sans Serif"))
par(mar=c(2.5, 2.5, 1, 1) + 0.1)
plot(x=Ta0, y=Pa1/1000, type="l", #isopleth of G=0 for MinAct
     main="", xaxt="n", #yaxt="n",bty="n", #family="A", 
     xlim=c(30,60), ylim=c(0,6), xaxs="i",yaxs="i", #start x y axes from origin 
     xlab=expression('Air temperature '*italic(T)[a]*' (°C)'), ylab=expression('Vapor pressure  '*italic(e)[a]*' (kPa)'), 
     lwd=1, lty=1, pch=15, col="black", cex.main=1, cex.axis=0.6, cex.lab=0.7, las = 1, mgp=c(1.2,0.5,0), tck=-0.02)
axis(side=1, cex.axis=0.6, tck=-0.02, padj=-2.5)
lines(x=Ta0, y=Pa2/1000, lwd=1, lty=1, col="black", cex=1.) #isopleth of G=0 for LightAmb

points(x=Ta_m, y=Pa_m/1000, pch=16, lwd=1, col="black", cex=0.8)
points(x=Ta2_m, y=Pa2_m/1000, pch=17, lwd=1, col="black", cex=0.8)
#standard deviation of Pa when Ta is fixed at 36, 38, 40
for (n in 1:3) {
  lines(x=c(Ta_m[n],Ta_m[n]), y=c((Pa_m+Pa_sd)[n],(Pa_m-Pa_sd)[n])/1000, lwd=1, lty=1, col="black", cex=1.)
  lines(x=c(Ta2_m[n],Ta2_m[n]), y=c((Pa2_m+Pa2_sd)[n],(Pa2_m-Pa2_sd)[n])/1000, lwd=1, lty=1, col="black", cex=1.)
}
#standard deviation of Ta when Pa is fixed at 2.7, 2.1, 1.6 kPa
for (n in 4:6) {
  lines(y=c(Pa_m[n],Pa_m[n])/1000, x=c((Ta_m+Ta_sd)[n],(Ta_m-Ta_sd)[n]), lwd=1, lty=1, col="black", cex=1.)
  lines(y=c(Pa2_m[n],Pa2_m[n])/1000, x=c((Ta2_m+Ta2_sd)[n],(Ta2_m-Ta2_sd)[n]), lwd=1, lty=1, col="black", cex=1.)
}
text(34,4.75,expression(italic(G)*" = 0"), srt=325, cex=0.6)
text(34,3.72,expression(italic(G)*" = 0"), srt=325, cex=0.6)
text(55.5,0.83,"MinAct", srt=322, cex=0.6)
text(50.5,0.75,"LightAmb", srt=322, cex=0.6)

#add contour of Tw isopleths
contour(x=Ta0,y=Pa0,z=Tw, levels = seq(15, 40, by = 5), lty=2, vfont = c("sans serif", "bold"), #col="brown",
        add = TRUE, labcex = 0.6, method = "flattest") 

dev.off()



#======redraw with RH on the Y axis=======
#contour of TW with Ta and RH
RH0 <- seq(0, 100, 5)
Tw2 <- matrix(nrow=length(Ta0),ncol=length(RH0))
for (n in 1:length(Ta0)) {
  for (m in 1:length(RH0)) {
    esat_air0=Esat.slope(Ta0[n], formula="Alduchov_1996")$Esat
    Pa0 <- RH0[m]/100*esat_air0
    VPD <- esat_air0 - Pa0 #kPa
    Tw2[n,m] <- wetbulb.temp(Ta0[n], pa/1000, VPD, Esat.formula="Alduchov_1996") #degC
  }}

#jpeg(paste0(wd,"/Supplementary_Fig.9.jpg"), units = "cm", width = 9, height = 7.5, res = 300)
pdf(paste0(wd,"/Supplementary_Fig.9.pdf"), width = 9/2.54, height = 7.5/2.54 )

par(mar=c(2.5, 2.5, 1, 1) + 0.1)
plot(x=Ta0, y=RH1, type="l", main="", #xaxt="n",yaxt="n",bty="n", #family="A", 
     xlim=c(30,60), ylim=c(0,100), xaxs="i",yaxs="i", #start x y axes from origin 
     xlab=expression('Air temperature '*italic(T)[a]*' (°C)'), ylab=expression('Relative humidity (%)'), 
     lwd=1, lty=1, pch=15, col="black", cex.main=1, cex.axis=0.6, cex.lab=0.7, las = 1, mgp=c(1.5,0.5,0), tck=-0.02)
lines(x=Ta0, y=RH2, lwd=1, lty=1, col="black", cex=1.)

points(x=Ta_m, y=RH_m, pch=16, lwd=1, col="black", cex=0.8)
points(x=Ta2_m, y=RH2_m, pch=17, lwd=1, col="black", cex=0.8)
#standard deviation of Pa when Ta is fixed at 36, 38, 40
for (n in 1:3) {
  lines(x=c(Ta_m[n],Ta_m[n]), y=c((RH_m+RH_sd)[n],(RH_m-RH_sd)[n]), lwd=1, lty=1, col="black", cex=1.)
  lines(x=c(Ta2_m[n],Ta2_m[n]), y=c((RH2_m+RH2_sd)[n],(RH2_m-RH2_sd)[n]), lwd=1, lty=1, col="black", cex=1.)
}
#standard deviation of Ta when Pa is fixed at 2.7, 2.1, 1.6 kPa
for (n in 4:6) {
  lines(y=c(RH_m[n],RH_m[n]), x=c((Ta_m+Ta_sd)[n],(Ta_m-Ta_sd)[n]), lwd=1, lty=1, col="black", cex=1.)
  lines(y=c(RH2_m[n],RH2_m[n]), x=c((Ta2_m+Ta2_sd)[n],(Ta2_m-Ta2_sd)[n]), lwd=1, lty=1, col="black", cex=1.)
}
text(34.1,90,expression(italic(G)*" = 0"), srt=300-5, cex=0.6)
text(32.6,73,expression(italic(G)*" = 0"), srt=300-3, cex=0.6)
text(52.6,6.8,"MinAct", srt=336, cex=0.6)
text(47,9,"LightAmb", srt=330, cex=0.6)

#add contour of Tw isopleths
contour(x=Ta0,y=RH0,z=Tw2, levels = seq(15, 40, by = 5), lty=2, vfont = c("sans serif", "bold"), #col="brown",
        add = TRUE, labcex = 0.6, method = "flattest") 

dev.off()
